function [peakIndex, peakValue] = speaks_movmedian(y, medianSlideLen)
	
	if ~isvector(y), error('Input must be vector'); end
	if iscolumn(y), y = y'; end
	
	flagPlot = false;
	
	dy = diff(y);
	maxIndex = find(diff(sign(dy))==-2) + 1;
	
	noise = movmedian(y, medianSlideLen); % 滑窗中值作为各处本底噪声的估计
	ldrValue = y(maxIndex) - noise(maxIndex);
% 	ldrs = y - noise; % return ldrs
	% 局部极大与本底噪声的差, 称为本地动态范围 LDR
		
	% 依据本地动态范围的大小排序, 进一步筛选出前 (1-r)x100% 分位的数据点
	r = .95;
	[~, idxIdx] = maxk(ldrValue, ceil((1-r)*length(maxIndex)));
	
	% 查看数据的 LDR 的分布特性
	if flagPlot
		figure; hold on;
		colormap('lines');
		yyaxis left
		histogram(y-noise, 20, 'DisplayName', 'raw data ldr');
		histogram(ldrValue, 'DisplayName', 'local max ldr');
		ylabel('counts');
		yyaxis right
		histogram(ldrValue(idxIdx), 'DisplayName', ...
			sprintf('local max ldr top 5%%: %d', length(idxIdx)), 'Normalization', 'probability');
		ylabel('ratio');
		legend
	end
	
	% 这些数据的横坐标按照所属的峰分为若干簇, 使用聚类方法自适应的确定个数
	index4clustering = maxIndex(idxIdx);
	[clustereLabels, ~] = dbscan(index4clustering', ceil(medianSlideLen/10), 1);
	
	peakCount = max(clustereLabels);
	if peakCount <= 0 % 聚类不成功
		warning('Fail to cluster');
		peakIndex = index4clustering;
	else % 聚类成功
% 		% 每簇内, 取簇的横坐标均值作为结果.
% 		peakIndex = zeros(1, clusterCount);
% 		for i = 1:clusterCount
% 			peakIndex(i) = ceil(mean(clusteringIndex(clustereLabel==i)));
% 		end
		% 每簇内, 取纵坐标最大点的横坐标作为结果.
		peakIndex = zeros(1, peakCount);
		for peakIdx = 1:peakCount
			clusterIdx = index4clustering(clustereLabels==peakIdx); % 提取这个簇内点的原始坐标
			[~, maxClusterIdx] = max(y(clusterIdx)); % 根据这些坐标对应的幅值求最大值, 作为这个簇的最终结果
			peakIndex(peakIdx) = clusterIdx(maxClusterIdx); % 提取最终结果的原始坐标
		end
	end
	peakValue = y(peakIndex);
% 	peakLDR = y(peakIndex) - noise(peakIndex); % return peakLDR
	
	% 查看频谱和 LDR
	if flagPlot
		
		figure('Units', 'normalized', 'Position', [0 0 .6 .8]);
		tiledlayout(2,1);
		
		nexttile; hold on;
		plot(y, 'DisplayName', 'Input Data');
% 		plot(maxIndex, y(maxIndex), '.', 'DisplayName', '局部极大');
		plot(noise, 'k', 'DisplayName', 'Noise Estimate'); % 滑窗中值

		% 峰值展示方法: 直接描点, 或依次描点并加入图例, 或在点附近标注
	 	plot(peakIndex, peakValue, 'd', 'DisplayName', '信号峰值');
% 		for i = 1:peakCount
% 			plot(peakIndex(i), peakValue(i), 'd', 'DisplayName', ...
% 				sprintf('%d, %.1f@%d', i, peakValue(i), peakIndex(i)));
% 			text(peakIndex(i), peakValue(i), ...
% 				sprintf('\\leftarrow %d', i));
% 		end
% 		plot(peakIndex, peakValue, 'd', 'DisplayName', '信号峰值');
% 		for i = 1:peakCount, text(peakIndex(i), peakValue(i), sprintf('\\leftarrow %d: %d\n    %.1f', i, peakIndex(i), peakValue(i))); end
		legend; % ('Location', 'northeast') eastoutside
		
		nexttile; hold on;
		title(sprintf('滑窗宽度 %d, 频谱宽度 %d', medianSlideLen, length(y)));
		plot(maxIndex, ldrValue, '.', 'DisplayName', 'Local DR');
		plot([1 length(y)], [1 1].*ldrValue(idxIdx(end)), '-.', 'DisplayName', '95%');
		legend;
	end
	
end